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ABSTRACT 



This paper describes the first steps in the development of a new multi-dimensional time implicit code devoted to 
the study of hydrodynamical processes in stellar interiors. The code solves the hydrodynamical equations in spherical 
geometry and is based on the finite volume method. Radiation transport is taken into account within the diffusion 
approximation. Realistic equation of state and opacities are implemented, allowing study of a wide range of the problems 
characteristic of stellar interiors. We describe the numerical method and various standard tests performed to validate 
the method in detail. We present preliminary results devoted to describing stellar convection. We first performed a local 
simulation of convection in the surface layers of a A-type star model. This simulation tested the ability of the code to 
address stellar conditions and to validate our results, since they can be compared to similar previous simulations based 
on explicit codes. We then present a global simulation of turbulent convective motions in a cold giant envelope, covering 
80% in radius of the stellar structure. Although our implicit scheme is unconditionally stable, we show that in practice 
there is a limitation on the time step that prevents the flow moving over several cells during a time step. Nevertheless, in 
the cold giant model we reach a hydro CFL number of 100. We also show that we are able to address flows with a wide 
range of Mach numbers (10 -3 < M s < 0.5), which is impossible with an anelastic approach. Our first developments 
are meant to demonstrate that applying an implicit scheme to a stellar evolution context is perfectly thinkable and to 
provide useful guidelines for optimising the development of an implicit multi-dimensional hydrodynamical code. 

Key words. Hydrodynamics - Convection - Methods: numerical - Stars: interiors 
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1. Introduction 

With the advent of massively parallel computers and the 
development of advanced numerical techniques like adap- 
tive mesh refinement and massive parallelisation, astro- 
physical fields like cosmology, star formation, or stel- 
lar/galactic environment studies have recently made re- 
markable steps forward. Despite this revolution in computa- 
tional astrophysics, our physical understanding of stellar in- 
teriors and evolution still largely relies on one-dimensional 
calculations. Implicit ID stellar evolution codes were widely 
developed during the past century, gaining more and 
more in sophistication through improvements in the in- 
put physics and implementation of an increasing number of 
complex physical processes such as time-dependent turbu- 
lent convection, rotation, or MHD processes. Their descrip- 
tion, however, still mostly relies on simplified, phenomeno- 
logical approaches, with a predictive power hampered by 
the use of several free parameters. These approaches have 
now reached their limits for understanding stellar struc- 
ture and evolution. The development of multi-dimensional 
hydrodynamical simulations becomes crucial for progress 
in the field of stellar physics and for the enormous ob- 
servational efforts aimed at producing data of unprece- 
dented quality, as expected from the asteroseismological 
space projects COROT and Kepler. 

Several efforts have been devoted to the development 
of 2D and 3D tools for studying some aspects and spe- 



cific problems of stellar structure and evolution. Impressive 
progress and developments have been achieved in multi- 
dimensional hydrodynamics and MHD models of stellar at- 
mospheres and photosp h eres (see e.g. Frevtag et al.l 19961: 
[ Stein fc Nordlundl Il998t lAsplund et alj 12000: iBigot et all 
2006; iNordlund fc Steinl I2009T ) . In stellar interiors, many 
studies of convection, rotation, and magnetic fields are 
performed with anelastic hydrodynamic solvers, which fil- 
ter out sound waves and have the priceless advantage of 
not being restricted to the Courant-Friedrich-Lewy (CFL 
hereafter) time step limit, a severe li mitation for ex- 
plicit compressible hydrodynamic solvers (jBrun fe Toomrel 
120021 ; iGuzikl [20101 and references therein). Anelastic ap- 
proaches, however, are restricted to the study of flows char- 
acterised by very low Mach numbers and small thermody- 
namic fluctuations from the background (see d iscussions in 
lAlmgren et~aT1l2006t Ovleakin fc Arnettll2007aD . 

Other groups use 2D/3D explicit hydrodynamical sim- 
ulations to follow the properties on short snapshots of 
turbulent convection, mixing processes, or nuclear burn- 
ing in stellar interiors such as cores and burning shells 
of massive stars or during the core He flash in low- 



mass stars dDearborn et al.l 20061: Mea 


kin & Arnetd 2007bl 


lEggleton et al. 120081: iMocak et al.l 


20091: lArnett et al.l 



of explicit codes, simulating those processes on time scales 
that are significantly larger than the dynamical one and 
relevant to stellar evolution would be prohibitively expen- 
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sive. To overcome this problem, increasing efforts are now 
devoted to developing algo rithms for low Mach-number , 
compressible stellar flows (|Huieirat fc Thielemann! l2009t 
iNonaka et aD l2010l ). However, many stellar regions and 
evolutionary phases are characterised by low-to-moderate 
Mach-number flows, which would ideally be described by 
fully time-implicit solvers. 

We only know two multi-dimensional implicit codes ap- 
plied to stellar physics problems, namel y the 2D stellar 
evolution code ROTORC developed b v iDeupred (|f990h 
and the 2D code VULCAN of iLivnd f|l993h . The work 
of Dcupree is pioneering in the field of multi-dimensional 
implicit stellar evolution, with ROTORC applied to the 
study of stellar convection and rotation in stars of differ- 
ent masses. However, the numerical method involved in this 
code cannot compete with current, advanced numerical al- 
gorithms, because it is based on a non-conservative finite- 
difference scheme with the implicit solver as an extension 
to the 2D of the Henyey method. This method is widely 
used in ID stellar evolution calculations but cumbersome 
(if not impossible) to parallelise, thus it significantly re- 
stricts the spatial resolution for multi-dimensional calcula- 
tions. VULCAN is based on a two-step method consisting of 
an implicit Lagrangian step followed by an explicit remap- 
ping step on an arbitrary grid. The code has been essen- 
tially applied to s t udying core-col l apse s upernovae (see e.g. 
iLivne et al.l l2004t Burrows et al. 2007) and novae explo- 



sions (sec c.6 



g. lGlasner fc Livnelll995[lLivne fc Arnetdll995t 



lOlasner et all 119971: iLivnd Il999h . Although VULCAN of- 
fers the possibility of solving advection implicitly, to our 
knowledge in these works the authors performed multi- 
dimensional computations by solving advection explicitly 
and radiation implicitly. Thus, we cannot gauge the ca- 
pacity of this code to describe long-term stellar evolution 
problems, in general, a nd convective flo ws in various stellar 
interiors, in particular . iHui eiratl (2005a. b) develops a multi- 
dimensional implicit MHD solver, but the applications are 
oriented more towards high-energy processes (e. g. accretion 
on a b lack hole) . Finally, we menti on the work of lToth et al.l 
(|l998f) and iKeppens et al.l (|l999h , which describe the im- 
plementation and tests of various implicit methods for 
multi-dimensional computations of steady state and time- 
dependent problem in (magneto)hydrodynamic. 

In this context, we started to develop a multi- 
dimensional time implicit code, guided by the motivation of 
improving the description of stellar interiors during differ- 
ent phases of evolution. The holy grail is to develop a tool 
that allows the 3D description of a complete star during 
thermal, nuclear, or accretion time scales that would be rel- 
evant to various stellar evolution phases. Our primary inter- 
ests focus on the description of time-dependent convection 
in the envelope of pulsating stars like Cepheids, a long- 
standing problem in asteroseismology, or during the very 
early phases of stellar evolution (early pre-main sequence) 
following the hydrodynamical phase of a proto-star forma- 
tion. Improving our description of these very early stages 
of evolution is crucial to a more accurate analysis of cur- 
rent observations and, in particular, to derivation of a key 
property for understanding star forma tion, namely the ini - 
tial mass function (see e.g discussion in lBaraffe et al.| [2002'). 
For the above-mentioned problems, neither an anelastic nor 
a low-Mach approach is appropriate, since the convective 
flows are characterised by very low (M <C 0.01) to moder- 
ate (O.f < M < 1) Mach numbers from the deep interior 



to the stellar surface. Also, a wide range of other astro- 
physicaly interesting problems can be addressed with such 
a tool, providing a wealth of applications on a long-term 
future and a revival of the field of stellar physics. 

Our goal is exceedingly challenging, with a long way 
to go before achievement. This first paper describes our 
method and the first tests performed to describe convec- 
tion in stellar interiors. Because implicit schemes are much 
more expensive in terms of CPU than explicit ones, part of 
the astrophysical, computational community casts doubts 
on the advantage and applicability of the former to any stel- 
lar physics problem. Our first developments are meant to 
demonstrate that the use of an implicit scheme in a stellar 
evolution context is perfectly thinkable. For this purpose, 
we present two examples of stellar convection for which we 
derive the exact CFL numbers and discuss the advantages 
and limits of an implicit approach. Even though the im- 
plicit method used throughout the paper has no stability 
limit on the time step (see Sect. 13.2^21) . the time step cannot 
take arbitrarily high values for both technical and accuracy 
reasons. We show in this paper that both reasons are ac- 
tually related and that the time step cannot grow much 
larger than the shortest time scale of the flow for cross- 
ing a grid cell. Technically this implies that the solution 
change is moderate between two time steps, insuring a safe 
convergence in the non-linear solver. In addition, this is 
also a good accuracy criterion, since advection over more 
than one grid cells is subject to strong numerical damping. 
Our first numerical experiments are meant to provide use- 
ful guidelines for developing an implicit multi-dimensional 
hydrodynamical code. 

The paper is organised as follows. We describe in Sect. 
[5] our physical model and in Sect. [3] our numerical method. 
We present a test case in Sect. |4] and results for stellar 
models in Sect. [5j In Sect. [51 we conclude the paper and 
discuss future numerical advancements that are planned to 
optimise the code. 



2. Equations 

We solve the equations describing the evolution of density, 
momentum, and internal energy, taking into account exter- 
nal gravity and radiative diffusion: 



dl P 



dl pU 



-V.(pu) 

-V.(peu) - PV.u - V.F r + —-, 
- V.(pH -r)- VP + pg, 



(1) 
(2) 
(3) 



where p is the density, e the specific internal energy, u the 
velocity, P the gas pressure, F r the radiative flux (see be- 
low), g the gravitational acceleration, and f the viscous 
stress tensor, whose components are given by 



Ipv 



(4) 



where is the strain rate tensor and v the viscosity co- 
efficient. The expected value of the molecular viscosity in 
stellar interiors implies large Reynolds numbers character- 
ising the flow (up to fO 12 ). It is therefore impossible to 
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model all scales of the flow, from the stellar scale down to 
the dissipation scale, on current generation of computers. 
As a result, any numerical simulations of stellar interiors 
should be interpreted in the large eddy simulation (LES) 
paradigm, in which only large-scale motions are resolved 
on the grid. In the standard LES, the effect of subgrid- 
scale motions (i.e. turbulence) is taken into account by in- 
troducing a so-called subgrid scale (SGS) model. A clas- 
si cal example is the ef fective viscosity coefficient proposed 
in Smagorinsky (1963). On the other hand, the MILES ap- 
proach (for Monotone In tegrated Large Eddy simulation, 
see e.g. iBoris et al.lll992f ) relies solely on the numerical vis- 
cosity of the scheme, which results from discretisation, to 
mimic the subgrid scale dissipation. In our case, a stan- 
dard LES approach should rely on an SGS model relevant 
for compressible flow in a stratified medium, which to our 
knowledge has not been designed yet. It is far from certain 
that the use of standard SGS models, mostly relevant for 
incompressible and homogenous flow, will improve our re- 
sults. Therefore, we choose to follow the MILES strategy 
and all the computations of stellar convection presented in 
this paper are based on the advection-diffusion equations 
without explicit viscosity. 

We solve the internal energy density equation. Some 
of the popular numerical codes that have been used to 
study stellar convection solve the entropy equation (e.g. 
ASH code, Pencil code), which is equivalent to the internal 
energy equation. A better approach to energy conservation 
would be to solve the total energy equation, as for instance 
in the COBOLD code (se e Freytag and collaborato rs) and 
in the PROMPI code fsee iMeakin fc ArnettlfeOOTbl ). In the 
future wc plan to implement a total energy equation solver, 
which might be more appropriate for problems like stellar 
pulsations. 

For stellar calculations, we treat radiation transport in 
the diffusion limit approximation. This approximation is 
suitable for an optically thick region but becomes inaccu- 
rate in the photosphere and optically thin regions. In this 
framework, the radiative flux F r writes as 



F r = ~k lad VT, 
where the photon conductivity fc ra d is given by 

4acT 3 



^rad 



3tzp 



(5) 



(0) 



with k the Rosseland opacit y of the gas, which is in terpo- 
lated from the OPAL tables (jlelesias fc Rogerslll99fih . 
The equations are closed with the equation of state: 



P = P(p,e), T = T(p,e) 



(7) 



Our equation of state is tabulated and includes the par- 
tial ionisation of several chemical elements (determined by 
the Saha-Bolztmann equations) from hydrogen to silicium 
and takes Coulombian pressure into account. The chemi- 
cal abundances are set to the solar values. For details see 
lAubert et al.l (|l996h . We consider here stellar envelopes and 
therefore we do not take nuclear reactions into account, 
which concern the central region of the star. 

For the geometry, we adopt spherical coordinates. We 
assume azimuthal symmetry so the only independent coor- 
dinates are r (the radius) and (the colatitude). In gen- 
eral, the coordinates (r, 9) span a domain [n n , r out ] x 2 ]. 



We use periodic conditions in the tangential direction. For 
boundary conditions at the top and bottom of the domain 
we impose 



1. non-penetrative conditions: u r — at r = ri n ,r ou t; 

2. stress-free conditions: -§p(^f-^j at r = rj n ,r out ; 

3. constant energy flux F* at the bottom of the domain 
r = r in ; 

4. an energy flux at the top of the domain given by 



(8) 



where T out is the temperature of the cells in the top layer of 
the domain. Finally, we would like to stress that we adopted 
2D geometry to simplify the development of the code, but 
extension to 3D is planned in the future. 



3. Numerical method 

We present below our numerical method, describing sepa- 
rately the spatial (Sect. 13. ip and temporal (Sect. I3.2[) dis- 
cretisations. 



3.1. Spatial discretisation 

3.1.1. One-dimensional advection scheme 

The equations are discretised on a staggered grid, using 
a finite volume approach. Cell interfaces location are de- 
noted by Xi,i — X...N X and cell centres by 2^+1/2, « = 
1 .. .N x — 1. In the staggered grid approach, scalar quan- 
tities (e.g. density, internal energy, pressure, temperature, 
etc) are located at cell centres, while velocity components 
are located at cell boundaries. 

In the finite volume approach, the physical equations 
are integrated over a control volume to yield the evolu- 
tion law for the physical quantities (mass, internal energy, 
momentum) in each cell. Fluxes are computed at cell inter- 
faces, yielding a conservative scheme with respect to advec- 
tion. For scalar quantities (density and internal energy) at 
x i+ i/2, the control volume is the cell [xi, li+i] with volume 
V i+1 /2 and surfaces Si, Si+i. For the velocity component 
at Xi, the control volume is the cell [^1-1/2)^1+1/2]; with 



volume Vi and surfaces Si-1/2, S^+i/a- 

We use the Van Leer method to co mpute the fluxes (see 
van Leerll977t lSt one Sz Normanll992|) . To compute the flux 
at the control volume interfaces, it is necessary to perform 
interpolation since not all quantities are available at the 
considered interface. We define "barred" quantities (i.e p, u) 
as quantities that are interpolated at the interface by tak- 
ing a simple average between the values available on each 
side of the interface, for instance Uj+1/2 = (itj + Uj+i)/2. 
For a stable and TVD advection, it is necessary to intro- 
duce some limiting process in the interpolation scheme of 
advected quantities. Here we proceed with an upwind lim- 
ited interpolation similar to the MUSCL method. In each 
cell a linear reconstruction of the solution is performed. For 
instance, we describe below the procedure for cell i + 1/2 
and cell centred variable q. We first compute the limited 
slope Aj_|_i/2<? of the linear reconstruction in the cell by 
using the values of the neighbouring cells with a limiting 
process: 
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Ai+iq = ?i + 3/2 - Qi+1/2 
= Qi+1/2 ~ Qi-1/2 



A;9 



(9) 
(10) 

(11) 



where <fi is the slope limiter. In the following we use the Van 
Leer limiter: 



2r 



if r > 



1 + r 
— 0, otherwise. 



(12) 
(13) 



With such a reconstruction scheme, we reconstruct at 



each interface right and left values q R and q L ' 



Qi+1/2 
Qi-l/2 



A i+ i/ 2 <? 

A l _ 1/2 g. 



(14) 
(15) 



Depending on the sign of the velocity, one of these quan- 
tities is used to compute the flux at the interface. For in- 
stance, at interface Xj one has 



<7* = q 



L if Ui > 



= qf- if Ui < 0. 



(16) 
(17) 



This defines the "tilded" quantities in the equations shown 
below. For the interfaces located at i + 1/2, the same re- 
construction is performed and the interface value is chosen 
depending on the sign of u i+ i/2- 

The pressure gradient at the interface is computed with 
a second-order, centred difference. The overall scheme is 
expected to be second order in space in the smooth regions 
of the flow. We present in appendix [X] an order study for 
the linear advection problem. 

After integration over the corresponding control volume 
and discretisation, the ID Euler equations become 



^A+l/2 e i+l/2^i+l/2 



dt 



- {Si+ipi+iu i+ i - SipiUi) (18) 

- (Si+ipe l+1 u. l+ i - S l pe i u t ) 

- Pi+yi{S%+iu i+ \ - SiU t ) (19) 

PiUiVi = - (S i+ i/2pUi + i/ 2 Ui + l/2 
~ Si-l/2P u i-l/ 2 Ui-l/ 2 ) 



Pi+1/2 — Pi 



Ax 



t ' 1/2 -V l . (20) 



3.1.2. Gravity 



Since we are solving the internal energy equation, gravity 
only enters the momentum equation. We assume that grav- 
ity remains constant and spherically symmetric so that the 
gravitational acceleration at radius T{ is given by 



GAL 



9i 



(21) 



j+1- 



j+l/2 - 





1 
1 

1 

1 






— 1 — 

1 p,e 
Ur 






1 
1 
1 





Fig. 1. Staggered mesh in 2D 



where G is the gravitational constant and Mi the mass 
contained inside the radius rj. Our stellar models do not 
extend to the centre of the star so Mi is computed assuming 
M\ = Mcoro, where M corc is determined from the original 
stellar model. 

Under the assumption of spherical symmetry, the grav- 
itational force enters the radial momentum equation only 
through the additional source term: 



PigiVi. 



(22) 



The assumption of spherical symmetry could become in- 
valid if the star developed significant non-radial oscilla- 
tions. The convective flows, which are not symmetric, do 
not alter the mass distribution symmetry significantly. 
Implementation of a full solver for the Poisson equation 
is planned for the future. 

3.1.3. Radiative diffusion 

The radiative flux through the cell interface is computed 
with a second-order central difference 



where 



1 



F r 



ac 1 1 t+l/2 1 i-l/2 



3 p^ 



1 



Ax 



1 



Pi+l/2 K i+l/2 Pi-l/2 K i-l/2 



(23) 



(24) 



3.1.4. Extension to 2D and spherical geometry 

Extension to 2D is done in the simplest possible way. Fluxes 
are computed along both directions following the ID algo- 
rithm sketched above and added together in the update step 
(see below). Figure [T] shows the location of all quantities on 
our staggered 2D mesh. We avoid dimensional splitting (e.g. 
Strang splitting) that would increase the cost of the solver 
since several (implicit) substeps would have to be computed 
to advance the solution over one time step. It is clear in this 
case that the accuracy of the method as deduced from ID 
analysis does not extend to multi-dimensional computation. 
One solution to improve the accuracy of the code would be 
to i mplement a genuinely mult i-dimensional method (see 
e.g. IColellalfl990t iLeVequell 19971) . Another possibility would 
be to increase the spatial order of the ID reconstruction 



4 



M. Viallet et al.: Implicit calculations of stellar interiors 



scheme to obtain a better accuracy in multi-dimensional 
computations. This is left for future developments. 

We describe below how spherical geometry is taken into 
account in our spatial discretisation: in our finite volume 
approach, the cell (i + 1/2, j + 1/2) spans the parameter 
space domain [rj,rj + i] x [9j,6j + i]. The geometrical factors 
in the 2D version of cell-centred equations reads as 



For the hydrodynamic equations, which are hyperbolic 
equations, the CFL condition states that the time step 
should not be longer than the crossing time of the fastest 
wave over a cell. In this case the maximum stable time step 
is equal to 



At 



Ax 



hydro 



(30) 



2ir 

Vi+i/2,3+1/2 = y Of+i - »i)(cos0j - cos0,-+i) (25) 



S ij+i/2 = 27rr?(cos0,--cos0,-+i) 
s i+i/2,j = 7r sin0 i (rf +1 - rf ). 



(26) 
(27) 



These geometrical factors account for the geometrical 
effect associated with spherical coordinates. There are sim- 
ilar expression for the radial momentum and the tangential 
momentum equations (see appendix[B|); however, due to the 
vector nature of the velocity, there are additional geometric 
source terms that appear in the momentum equations. In 
the radial momentum equation, the geometric source term 
reads as 



ogcom _ -Tr 



Pi, ] + l/2{u 9 ij+1/2 ) 2 



(28) 



and the geometric term in the tangential momentum is 



ogcom _ r> Pi+l/2,j u i+l/2, 3 U t+l/2J , , 

i+l/2,j - - V i+l/2,j ; ■ (^9 J 

' r i+l/2 

The full semi-discretised formulae are detailed in the ap- 
pendix [Bj 

Spherical coordinates have severe drawbacks in the dis- 
cretisation of a sphere. First, since r = is a singular point 
of the coordinates, the cell size reduces to zero toward the 
centre. They are therefore inappropriate for modelling the 
central regions of the star. Also, the cell size increases to- 
ward the surface, wher eas finer resolut i on is n eeded in these 
regions (see Sect. 15.21 1. ICalhoun et al.l (|2008l) shows how to 
design structured grids that overcome these limitations. 



where \u\ + c s is the velocity of the fastest wave (u is the 
velocity of the background flow). 

When time-dependent diffusion processes are consid- 
ered, for instance in our case radiative diffusion, the cor- 
responding CFL time step is related to the typical time 
scale for diffusion in a cell: 



At 



rad 



Ax 2 

X 



(31) 



where x = &rad/(pc p ) is the radiative diffusivity coefficient, 
and c p the specific heat at constant pressure. For advection- 
diffusion problem, which are of mixed nature, the CFL time 
step is taken to be the minimum values of equations Q30p 
and ([51). 

We also define here a useful time step that we call the 
advective time scale. It is defined by 



At 



Ax 



adv 



(32) 



Its definition is inspired by the hydro time step (|30[) where 
we have set the velocity of the wave to be simply u instead 
of u + c s . This time step characterises the time needed for 
the flow to cross one cell. We stress that, for a code solving 
the standard compressible hydrodynamic equations with an 
explicit method, this time step is not a stability limit. We 
introduce this time scale since, as shown later, it corre- 
sponds to a natural limit for the time step; however, for 
anelastic codes, this time step represents the stability limit, 
since sound waves are filtered out in this approach. 

With the definition of these three time steps, we can de- 
fine the corresponding CFL numbers for a given numerical 
time step At: 



3.2. Temporal method 
3.2.1. CFL time steps 

When an explicit method is used to integrate a discrete 
set of equations in time, there is a constraint on the time 
step that results from the CFL condition. The CFL con- 
dition states that the physical domain of dependence of 
the solution should be included in the numerical domain 
of dependence, otherwise the numerical scheme cannot be 
convergent. For an explicit scheme, this yields a constraint 
on the time step that has to be smaller than the CFL time 
step to perform a stable computation (independently of any 
accuracy consideration) . 

Determining the stability limit for a given explicit 
scheme is not always trivial, if not impossible, especially 
in the non-linear and multi-dimensional cases. See, for in- 
stance, iHirschl (|1990f ) for a careful derivation of the CFL 
time step for popular schemes. The CFL time steps that 
we define in this section are based on very simple argu- 
ments and should only be considered as a guideline for the 
maximum stable time step. 



CFLhydro — -rr ; CFL rac i — — — ; CFL a d v — -jr. • 

^^hydro ^^rad ^^adv 

(33) 

By adopting an implicit integration strategy, our goal 
is to release the stability limit on the time step. Therefore, 
to appreciate the performance of our solver, we carefully 
monitor the values of these three CFL numbers. 

In 2D, the time steps defined above are computed in- 
dependently in both directions and the CFL time step is 
defined as the minimum of these values. 

3.2.2. Implicit strategy 

After the spatial discretisation method described in the pre- 
vious section, one obtains a system of coupled ODEs (see 
Appendix [B]) : 

^jf- = R(U ktl ), (34) 

where Ut,i is a vector gathering all the quantities appearing 
below time derivatives (i.e. p, pe, pu r , and pug, multiplied 
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by the appropriate volumes), and R(Uk,i) contains all fluxes 
and sources s: 

R(Uk,i) = — (SX+i,ji J fc+i,j - Sk,l F k,l 

+ s l,i+i F k.i+i ~ S k,i F k,i) ( 35 ) 

+Sk,lVk,l. 

To obtain an update formula we now apply an implicit 
temporal method to integrate (Eq. l34|) in time: 

U lT - U iJ = Ai(/3iT +1 + (1 - P)R n ) (36) 

where R n (resp. R n+1 ) is the r.h.s. of equation ([34]) eval- 
uated at time n (resp. n + 1). The first-order accurate 
Backward Euler method corresponds to (3 = 1. The second- 
order accurate Crank- Nicholson corresponds to f3 = 1/2. 
Both the Backward Euler and the Crank-Nicholson meth- 
ods are A-stable, which implies that they are uncondition- 
ally stabl^B For stiff equations, a more r elevant propert y 
is the L-stability of the scheme (see e.g iLeVeaud [2007). 
The Backward Euler method is L-stable but the Crank- 
Nicholson method is not. This implies that the latter can 
have difficulties with stiff problems and produce spurious 
oscillations. Indeed, when computing stellar models with 
high values of the radiative CFL (see the A-type star model 
in Sect. 15. we have found that it is more stable to use 
a value of (3 that is slightly higher than 1/2, for instance 
(3 = 0.55. In the future, we plan to implement the second- 
order backward differentiation formula, which is a second- 
order implicit multistep method having the A & L-stability 
property. 

Equation (|36[) describes a system of non-linear equations 
that defines the solution at the new time step, U n+1 . We use 
p, e, u r , and Ut as independent variables, so the quantities 
that appear in the time derivatives in our equations are 
thus composite quantities of the primitive variables. 

3.2.3. Newton-Raphson procedure 

To solve for the new time step p n+1 , e n+1 , u™ +1 , iig + , one 
has to solve the set of non-linear equations F(U n+1 ) = 0, 
where the residual F is defined as 

F(U n+1 ) = U n+1 -U n - M(f3R n+1 + (1 - j3)R n ). (37) 

Here, we proceed with a Newton-Raphson procedure. The 
procedure to find U n+1 is initialised by taking the current 
solution as an initial guess for the new solution: 

C/ (0) = U n . (38) 

At each Newton-Raphson iteration we solve for the correc- 
tions SU^: 

J x 5U {k) = -F(U (k) ), (39) 

where J is the Jacobian matrix of the residual vector F. 
The current iteration is then updated with 

jjik+i) =u {k) + x8UW, (40) 

1 Not all implicit methods are A-stable; e.g. the Adams- 
Moulton methods are implicit but only conditionally stable. 



where < A < 1 is introduced to prevent divergence of 
the procedure. The classical Newton-Raphson update corre- 
sponds to A = 1. A straightforward line-searching algorithm 
is applied to find the value of A that ensures a decrease in 
the residual. However, we found that when solving the non- 
linear radiative diffusion term at very high radiative CFL 
numbers (~ 10 6 ), the first Newton-Raphson iteration al- 
ways increases the residual. Different tests have shown that 
this seems unavoidable, and since the system would even- 
tually converge within a few iterations, we always enforce 
A = 1 during the first iteration in this particular case. 

We stop the iteration procedure when the relative cor- 
rection drops below a certain level: 



where Uq are typical values of the variables. For the den- 
sity and internal energy we simply choose the correspond- 
ing values at iteration k, but for the velocities we choose 
max(u, c s ), where u is the value of the velocity at iteration 
k, and c s is the sound speed. In our computations we set 
e = 10- 6 . 

3.2.4. Jacobian matrix computation 

The Jacobian matrix elements are the partial derivatives 
of the numerical scheme equations ([37|) with respect to all 
numerical variables. We compute the partial derivatives by 
finite differencing. The location of non-zero elements of the 
Jacobian depends exactly on the stencil of the numerical 
scheme (i.e. which variables contribute to a given equa- 
tion). It is therefore easy to determine the non-zero struc- 
ture (also known as the sparsity pattern) of the Jacobian 
matrix. Any Jacobian computation technique should then 
take advantage of the non-zero pattern to compute the 
Jacobian matrix elements. We use here the Co loured Finite 
Differencing technique (CFD hereafter; see ICurtis et al.l 
Il974t iGebremedhin et afll2005h . which minimises the num- 
ber of function evaluation needed to compute the Jacobian 
matrix. The CFD algorithm uses the non-zero pattern of 
the Jacobian to group independent columns (i.e. variables) 
of the Jacobian matrix, i.e. those that do not share a non- 
zero element on the same line, into "colours" . This yields a 
"column compressed" representation of the Jacobian. The 
number of colours n g , which is also the number of columns 
of the compressed Jacobian matrix, is roughly equal to the 
maximum number of non-zero elements found in a row. 
This number depends only on the stencil of the scheme and 
is (ideally) independent of the dimension of the Jacobian 
matrix. To compute the matrix elements, the CFD algo- 
rithm perturbs all variables associated with the same colour 
and recomputes the non-linear equations. The matrix ele- 
ments are then obtained by a straightforward finite differ- 
ence against a reference value of the equations. The num- 
ber of non-linear function evaluations is therefore equal to 
N{n g + 1), where N is the size of the matrix. As n g is in- 
dependent from N, the CFD algorithm is expected to scale 
as N. 

In our code, we use the CFD algorithm that is i mple- 
mented in the Trilinos library (see lHeroux et al.ll2003lh The 
graph-co louring algor i thm is the "greedy" algorithm pro- 
posed in ICurtis etail (|1974h . We find that for our matri- 
ces this algorithm typically produces a number of colour 
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Table 1. Run parameters for the entropy bubble test. 




Degrees of freedom (N) 

Fig. 2. Time spent in computing the Jacobian matrix 
(crosses) and in solving one linear system with MUMPS 
(dots) for different matrix sizes on a Pentium III Xeon 
quad-core CPU at 3 GHz. 

uq = 51 ± 1 for matrix sizes ranging from 2304 (20 2 do- 
main) to 652 864 (400 2 domain). Therefore, the number of 
colours is roughly independent of the matrix size, and the 
algorithm scales as N. This is illustrated in Fig. [5J where 
we plot the time needed to compute the Jacobian matrix 
for different spatial resolutions. 

3.2.5. Linear solver 

The Newton-Raphson procedure requires the solution of 
the linear equation (|3U)) . Equation (j3"9")l is a sparse linear 
problem that can be solved by either a direct method (i.e. 
LU decomposition) or an iterative method. Iterative solvers 
(e.g. GMRES, BiCGStab) perform usually faster than di- 
rect methods but they rely heavily on preconditioning to 
improve (or even obtain) convergence. On the other hand, 
direct solvers are more accurate and robust but they be- 
come very expensive as the size of the linear system in- 
creases. For convenience, and since our 2D approach al- 
lows it, we decided to use a direct solver for developing the 
code. We use in our code the MUMPS solver (MUlt i fronta l 
Massively Parallel Solver, see lAmestov et al~ll2001l [2006). 
MUMPS is able to handle general non-symmetric matri- 
ces; i.e., it is not restricted to tridiagonal or pentadiagonal 
matrices. MUMPS achieve robustness by detecting null piv- 
ots during the elimination phase and accuracy by applying 
a few (usually 2 or 3) steps of iterative refinement in the 
postprocessing phase. 

We use the advantage of knowing the non-zero structure 
of the Jacobian matrix to perform the symbolic factorisa- 
tion of the linear system only one time, at the beginning 
of a simulation. Later on, only the numerical factorisation 
is performed with a given r.h.s. vector to obtain the solu- 
tion of the linear system. The cost of the direct solver is 
illustrated in Fig. [2] where the time needed to solve one 
linear system is shown for various spatial resolutions. The 
linear solver dominates the time needed to compute the 
Jacobian matrix. With a simple linear regression in the 
log-log space, we find that in our case (i.e. for the Jacobian 
matrices that result from our temporal and spatial discreti- 
sation) the time needed to solve one linear system scales as 
TV 1 - 38 . For 2D runs with resolution up to 500 2 , the num- 
ber of variables is under one million. Linear systems are in 



Run 


Final time 


Resolution 


time step 


CFL 


1 


200 


50 2 


0.04 


1 


2 


200 


50 2 


0.5 


12.5 



this case still tractable with direct methods. The memory 
requirement for the LU decomposition ranges from 10 Mb 
for 2304 degrees of freedom to 9.4 Gb for 652 864 degrees of 
freedom. It is clear that larger systems (e.g. 3D problems) 
will require too much memory and computational time to 
be tractable. Therefore, we plan in the future to implement 
iterative solvers that are less expensive both in terms of 
memory and CPU time. 

3.2.6. Time step control 

In an implicit method, one has to introduce a time step 
control strategy, usually based on empirical grounds. We 
found that the following three strategies are useful, in the 
order of the most to the least basic: 

Convergence-based: If the Newton-Raphson procedure 
converges within five iterations, increase the time step 
by a factor of 1.5. If NR converges in more than ten 
iterations, decrease the time step by a factor of 2. 

Relative change- based: Compute relative change in the 
solution between time step n and n + 1. If it is lower 
than 10 %, increase the time step by a factor of 1.5. If 
it is larger than 20 %, decrease time step by a factor of 
2. 

Advection-based: Set the time step so that CFL a d v = 1. 

The last strategy is physically motivated since it en- 
sures that the fluid does not move across more than one 
cell width during one time step. For the preliminary results 
presented in this paper, we have used the first, least strin- 
gent, and least physical, strategy. In the future we plan to 
compare these three strategies, both in terms of accuracy 
and computational time. Finally, in some cases, we also 
provide maximum/minimum values for the CFL number or 
the time step. 

4. Test: oscillations of an entropy bubble 

Basic tests (see appendix [5} have been successfully done 
for code validation during its first stages of development, 
so we focus here on a test that is more relevant to the phys- 
ical processes we would like to inve stigate. It was taken 
from iDintrans &; Brandenbur d I2004J (DB hereafter) . We 
considered an isothermal stratified atmosphere with con- 
stant gravity perturbed by an entropy bubble. We com- 
puted the oscillations of this bubble around its equilib- 
rium position, and we analysed the spectrum of internal 
gravity and sound waves driven by the bubble. The ini- 
tial setup, parameters, and units normalisation are similar 
to DB and are therefore not reproduced here. The Brunt- 
Vaisala frequency of the atmosphere is N ~ 0.82. For this 
test, Eq. (1-3) were solved with the viscous terms included. 
The kinematic viscosity coefficient v was set to 5 x 10~ 4 
in all runs. The thermal conductivity is set to zero, as we 
were interested in the adiabatic evolution of the system. 
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At = 0.04 At = 0.5 
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Fig. 3. Accoustic modes in the (z,u>) plane for At — 0.04 (left) and At = 0.5 (right). 



We used a Cartesian version of the codffl and considered 
a domain (x,z) G [—0.5,0.5] x [0,1]. Periodic boundary 
conditions were used in the horizontal directions and non- 
penetrative, stress-free conditions were used at the bottom 
and top boundaries. All the run parameters are summarised 
in Table Q] For each run we set the time step to a constant 
value. 

As in DB, we describe each wave mode by two integers 
I, n > 0, which are the number of nodes in the horizontal 
and vertical directions, respectively. In this test we looked 
at two types of waves: 

1. Vertical sound waves that have frequencies lj — (n+l)7r. 
These waves have I = since they have no horizontal 
structure. 

2. Internal gravity waves that have frequencies lu < N and 
are characterised by / ^ (internal gravity modes can- 
not propagate in the vertical direction). 

See DB for the derivation of these properties. These 
two kinds of waves occupy well-separated regions of the 
frequency spectrum. Vertical sound waves have vertical pe- 
riods P < 2 and will therefore be more sensitive to the time 
step. Internal gravity modes have longer periods P > 7.3 
and will therefore be less sensitive to the time step. This 
test thus provides a good analysis of the influence of the 
time step on the accuracy of the numerical results. 

We first analysed the vertical sound wave spectrum by 
using method 2 described in DB, for runs 1 & 2 (cf. Table 
[1]). In this method, the wave spectrum is represented in the 
(z,lu) plane. For any given time step At, the Nyquist the- 
orem states that no periodic signal with frequency higher 
than ojjv = tt / At can be resolved. Also the temporal Fourier 
transform has a frequency resolution that is Aw = 2ir/T, 
where T is the duration of the simulation. The results are 
shown in Fig. El for two different time steps. For At — 0.04, 
we recognise the signature of the n = 0,1,2 modes. The 
n = 3 mode is below the range of the isocontour levels. 
Furthermore, all the modes are located at the correct fre- 
quencies (within Aw ~ 0.125). For At = 0.5, one has 
wjv = 2-7T thus normally allowing only for modes n = 
and n = 1. In the map, one can recognise the signature of 
the n = 0, 1, 2 modes. None of them is located at the correct 



2 Cartesian coordinates are used in this particular case be- 
cause of analytical solutions in this geometry. 



frequency, and it is surprising to find the n — 2 mode whose 
eigenfrequency is above the Nyquist frequency. In this case 
it is clear that the time step is too large to accurately re- 
solve the sound waves. To analyse the gravity modes, we 
used the method developed in DB: we projected the veloc- 
ity field on the anelastic eigenvectors. With this method 
we obtained the individual mode amplitudes q„ which are 
the coefficients in the eigenvector expansion. With such a 
decomposition it is straightforward to obtain the mode fre- 
quency. We checked that, for both runs, the gravity mode 
frequencies match the analytical predictions to the percent 
level. 

This test highlights the importance of the choice of the 
time step to correctly describe the physical process of inter- 
est. In the present case, sound-wave amplitudes are much 
lower than the amplitudes of the gravity modes (the initial 
setup is almost in hydrostatic equilibrium), and therefore 
they do not play an important role in this problem. One 
can then safely use a higher time step (e.g. At = 0.5) to 
study the internal gravity modes, as shown by the excellent 
agreement with predicted frequencies. Obviously, one could 
not do that if the initial setup was intended to be far from 
hydrostatic equilibrium, since sound waves would at least 
be as important as internal gravity modes. 



5. Stellar models 

We now test our code on realistic stellar conditions. 



5.1. Local simulation: A-type star 

We present in this section results of convection in A-type 
stars (see model I in Table [2]), since comparison with pre- 
vious models of such st ars, calculated by d i fferent groups 
using explicit codes (see lSofia fc ChaiJll984 iFrevtag et al.l 
1996), are possible. These types of stars are weakly con- 
vective and possess two convectively unstable zone near 
their surface: around T ~ 40 000 K due to Hell ionisation 
and near T ~ 10 000 K due to Hel/H ionisation. Usually 
the Hel/H surface convective zone is strong, but very nar- 
row, whereas the Hell convective zone is broader but much 
weaker. To model these convective regions it is therefore 
sufficient to focus on the surface layers of the star by per- 
forming local simulations. 
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Table 2. Stellar parameters of the models considered in this paper. 



Model 


Mass (M Q ) 


R (R@) 


T cff (K) log( ff ) 


Luminosity (Lg) 


Comments 


I 


2 


2.1 


8 500 4.10 


20 


A-type star 


II 


5 


58.8 


4 500 1.6 


1,258 


Cold giant 



Table 3. Numerical parameters for the A-star run. 



Resolution 


Domain 


TKH 


'-"hydro 


^^rad 


Final time 


Time steps 


Wall Time 


Mem. 


214 x 512 


[0.974, 1.] x [88.2°, 91.8°] 


1.15 d 


1.84 S 


3 x lO -2 s 


3.5 d 


4032 


2 w 5 d 


4.8 Gb 



t= 3.873e + 00 days 




-0.03 -0.02 -0.01 0.00 0.01 0.02 0.03 

Fig. 5. Snapshot of our simulation of an A-type star. The axes are Cartesian coordinates normalised by the stellar radius 
and with the origin at the star centre. 
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Fig. 4. Evolution of the total kinetic energy (top panel) and 
CFL numbers (lower panel, solid: radiative; dashed: hydro; 
dotted: advection) for the A-type star run. The vertical 
dashed lines show the time interval on which the averaged 
fluxes shown in Fig. [7] are computed. 



Our initial model starts at the photosphere, i.e. the lo- 
cation where t = 2/2, and T e s — 8 500 K. This particular 
choice inhibits the convection in the Hel/H ionisation zone, 
which is located at the photosphere, so we do not expect to 



model this convective region correctly. A better modelling 
of this region should include a part of the stellar atmo- 
sp here in the computati onal domain, as is done for instance 
by iFrevtag et al.l (|l996l) . Our initi al models are produced 
with a ID stellar structure code (see lBaraffe fc El Eidlll991|) 
where the mixing length theory (MLT hereafter) prescrip- 
tion for convection has been turned off. The resulting fully 
radiative models are very close to the MLT models (even 
near the surface) since convection is too weak, in terms of 
energy transport, to significantly modify the structure. The 
model presented here extends down to a depth where the 
temperature is T = 100 000 K, where the star is radiatively 
stable, so that the numerical domain includes a stable, ra- 
diative, buffer zone at the bottom of the domain. Since our 
simulation focusses on a very small fraction of the star, the 
geometry is Cartesian- like (but we use spherical geometry) . 
In terms of physical dimensions, the physical domain has 
an average width of 90 000 km and a height of 37 500 km. 

To construct the initial model for the multi-dimensional 
code, we simply copy the ID structure along the tangential 
direction. However, the number of tangential grid points 
has to be carefully chosen in order to keep an acceptable 
aspect ratio of the cells. The number of grid points in the 
tangential direction therefore depends on the number of 
grid points in the radial direction through the radial extend 
of the cells. 

We discuss below the results of our run. Table [3] sum- 
marises important parameters, tkh is the thermal time 
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Fig. 6. Radial profiles of the CFL time steps for the A- 
type star run (see text). From bottom to top: radiation 
(solid red), hydro (dash blue), and advection (dot green). 
The horizontal dash-dotted line indicates the time step of 
the snapshot under consideration. 
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Fig. 7. Radial profiles of averaged fluxes for the A-type star 
run (dashed line: enthalpy flux, dotted line: kinetic energy 
flux). The solid line is the convective flux as predicted by 
the mixing length theory with a mixing length parameter 
a = 1.375. The box in the upper left corner is a zoom of 
the surface layer of the star. 



scale defined as the ratio between the total internal energy 
in the computational domain divided by the luminosity of 
the model. The last column of Table [3] shows the mem- 
ory requirement for the LU decomposition. Calculations 
were done on an AMD "Istanbul" processor (6 cores at 
2.7 Ghz). Figure [4] shows the evolution as a function of 
time of the total kinetic energy in the domain (used to 
track the development and saturation of the convective in- 
stability) and the evolution of the different CFL numbers 
defined in Sect. 13.2.11 At the beginning, the model under- 
goes a relaxation toward hydrostatic equilibrium, and the 
CFL number rapidly reaches values > 10 5 . Eventually, the 
velocity field that develops as a result of the relaxation pro- 
cess triggers the convective instability, which corresponds 
to an increase in the kinetic energy to a significant level. As 
the first vortices develop, the time step decreases since the 
Newton-Raphson procedure would not converge otherwise. 
After ~ 2.3 days of simulated time, we reach a quasi-steady 
state that lasts until the end of the simulation at t ~ 3.5 
days. In this phase, the time step remains roughly constant 
with At = 18±2.6 s, corresponding to CFLhydro = 9.8±1.4 
and CFL ra( j = 400 ± 60 (the second value indicated corre- 
spond to the standard mean deviation). In this quasi steady 
state, convection is fully developed, and we observe a com- 
plex vortex dynamics involving vortex pairing/merging and 
the formation of downdrafts. These phenomena are charac- 
teristic of turbulent convection (see snapshot in Fig. [5]). 

Figure [6] shows the radial profiles of the different CFL 
time steps defined in Sect. 13.2.11 The CFL time steps are 
first computed in each cell of the domain, and the radial 
profiles shown in Fig. [6] are then deduced by taking the 
minimum value along the tangential direction for each ra- 
dius. We first note that the radiative CFL time step pro- 
vides the most stringent constraint. Solving the radiative 
diffusion with an explicit method would require a time step 
lower than Ai = 3xl0~ 2 s and would make any calcula- 
tion cumbersome. The hydro CFL time step increases with 
radius, so that the hydro limitation on the time step would 
come from the bottom of the numerical domain. This is ex- 



pected since we are using a uniform grid, and at the bottom 
of the envelope, the temperature, and hence the velocity of 
the sound waves, are the largest. The advective time scale 
shows the region where the flow velocities are the largest, 
i.e. near the surface. The time step of the actual simulation 
clearly corresponds to the lowest value of the advective time 
scale. If the time step was much longer than the advective 
time scale, the solution would change significantly since ed- 
dies would be moving across more than one grid cell. In 
this case, however, the Newton-Raphson procedure would 
have convergence problems since the initial guess (the solu- 
tion at time step n) would strongly depart from the correct 
solution. Since our time step strategy is based on limiting 
the number of iterations, the time step would be automat- 
ically decreased. Therefore, we do not expect our current 
implicit method to work with much larger time steps, and 
we identify the minimum value of the advective time scale 
as a natural (physical) limit for the time step. 

We now analyse the properties of the fully developed 
convective state. We computed the average of the radiative, 
mass, kinetic, and enthalpy fluxes: 

F rad = ({S r (-k^)) e )t 
ar 

F m = {{S r pu r )e)t 
F ki „ = ((S r u r ^pv 2 ) e ) t 

F on th = ((S r u r (pe + P))g-((pe + P))gF m ) t . 

These quantities have the dimension of a flux multiplied 
by a surface (i.e. a luminosity), in order to take geometrical 
effects introduced by the spherical coordinates into account. 
We however use the term "flux" in the following to make 
the discussion clearer. The total flux is defined as the sum 
of the enthalpy, radiative, and kinetic energy fluxes. 

Figure [7] shows the radial profile of the averaged fluxes 
defined above. The averaging was done on the interval 2.3 
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Fig. 8. Properties of the sampled grid for the giant model. Left panel: Pressure scale height to radial cell size ratio. Right 
panel: Cell aspect ratio when the radial mesh is extended in the tangential direction. 



days - 3 days (represents 1408 time steps) but we checked 
that, by taking smaller windows (but still large enough to 
have several hundreds of time step), we would obtain simi- 
lar results. Figure [7] shows the typical profile of the enthalpy 
flux in the Hell convective zone, as expected from general 
properties of convection: an inward/outward directed flux, 
with an outward part stronger than the inward flux. Also 
the kinetic energy flux is typical of what is expected from 
convection: downward, as a result of the acceleration of the 
downdrafts by gravity. The outward convective flux reaches 
a maximum amplitude of about 0.68 % of the stellar flux. 
Depending on the time interval on which the averaging was 
done, this maximum value ranges between 0.6 % and 0.77 
% of the stellar flux. This corresponds to the prediction of 
MLT for an a parameter in the range 1.35 — 1.4. As dis- 
cussed above, we do not expect to have realistic results for 
the H/Hel top convective zone since our outer boundary 
condition at the photosphere inhibits convection in this re- 
gion. However, despite the lack of resolution (only a few grid 
points resolve the ionisation region), we still find a non zero 
average convective flux with an amplitude of ~ 10~^F+. For 
the aforementioned value of a, tuned to match the Hell 
zone, the MLT predicts a convective flux in the Hel/H re- 
gion, which has a location and an amplitude consistent with 
our results. 

We now compa re o ur results w it h th e work of 
ISofia fc Char] f|1984T ) and iFrevtag et all (|1996t) . who have 
performed numerical models of a A-type star with the same 
parameters as ours (model A5 of Sofia & Chan paper, model 
AT85g41N3 of Freytag et al. paper). Sofia & Chen also con- 
sidered models extending to the photosphere, but they did 
not find the small convective flux that appears near the sur- 
face in our simulation, most certainly because of the lack of 
resolution in their comput ations (their resolut ion is twice 
lower than our). Since onlv lFrevtag et all (jl996l) resolve the 
Hel/H region properly, we restrict t he comparison to the 
Hell convective zone. For this region, ISofia fc Chan! (|1984[ ) 
found a convective flux between 0.15 and 0.2 % of the stel- 
lar flux (depending o n the aspect ratio of their model) and 
IFrevtag et al.l (|1996|) found a convective flux about 0.12% 
of the stellar flux. This is significantly lower than the am- 
plitude deduced from our simulations. This difference likely 
results from different input physics (metallicity, EOS, opac- 
ities) . 

Our results on A-type stars show that we could suc- 
cessfully describe 2D convective patterns in realistic stellar 



conditions with our fully implicit scheme. Thanks to this 
approach, the time step in our computation is not limited 
by sound waves or by radiative diffusion. However, we have 
shown that the time step is limited by the velocity of the 
fluid itself, which determines the rate of change in the solu- 
tion. The hydro CFL cannot be much larger than 10 when 
convective motions are fully developed. On the other hand, 
the radiative CFL number is very high (~ 400). This is ex- 
pected for the local models of "hot" stars that have very 
short radiative time scales in the low-density, photospheric 
region. The rather low value of the hydro CFL number 
and the high value of the radiative CFL number suggest 
that, for this type of local model, a mixed approach (ex- 
plicit hydrodynamic + implicit radiation) could allow sub- 
stantial speedup to be reached as compared to a fully ex- 
plicit/implicit approach. However, our goal with the present 
simulation of A-type stars is not to compete with previous 
simulations performed with explicit codes, which are doing 
an excellent job for this type of star with very tiny atmo- 
spheric and subphotospheric convective zones. Also a real- 
istic description of convection requires a 3D approach. The 
present study is meant to test the behaviour and the relia- 
bility of the implicit code to model convection under realis- 
tic conditions that have been already studied. Nevertheless, 
we stress that a fully implicit scheme remains very efficient 
at computing the initial development of convective insta- 
bility, since the time step can grow to high values as long 
as velocities remain low. 

5.2. Global simulation: cold giant star 

The next step is to compute models involving a significant 
fraction of the star. We consider here a model with pa- 
rameters corresponding to a cold giant star (model II in 
Table©. 

5.2.1. Initial model and grid design 

For this model, MLT predicts that about 50 % of the stel- 
lar envelope is convective. The initial model for the multi- 
dimensional simulations is based on a ID MLT model. 
Starting from a radiative model, as previously done for the 
A-type star, yields difficulties since a radiative structure 
departs strongly from the convective solution. The MLT 
model might also depart from the real convective state, 
but one can expect this difference to be less pronounced 
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Fig. 9. Properties of the uniform grid for the giant model with N r = 217. Left panel: Pressure scale height to radial cell 
size ratio. Right panel: Cell aspect ratio when the radial mesh is extended in the tangential direction. 



than with a fully radiative model. We model a significant 
fraction of the star (80 %) in order to have a large radia- 
tive zone below the convective zone. Our numerical domain 
therefore extends radially from 0.2i?* to i?*, where i?* is the 
radius of the star as defined by the location of the photo- 
sphere. This yields a pressure contrast of roughly six orders 
of magnitude (~ 14 pressure scale heights) and a density 
contrast of almost five orders of magnitude throughout the 
whole domain. Throughout the convective zone the pres- 
sure contrast is ~ 400 (~ 6 pressure scale heights) and the 
density contrast is ~ 50. Previous works that also consider 
an exte nded radial domain a re the "star- in-a-box" si mula- 
tions of iFrevtag et all ()2002f ) . IWoodward et al.l ([20031 ) . and 
iDobler et al.l ( 2006 ), where the full st ar is embedded in a 
Cartesian domain. TDobler et al.l (|2006f ) consider a very sim- 
plified stellar model, with a very low-density stratification 
(a factor of ~ 10). The two other works consider more re- 
alistic stellar models of super giant and AGB stars and 
obtain conve ction over multip l e pres sure scale heights. We 
also mention iBrun &: Palaciosl (|2009f ). which presents a nu- 
merical simulation of the inner ~ 50 % of a red giant star 
with the anelastic code ASH. Their density stratification in 
their convective zone (which occupies the whole numerical 
domain) is about two orders of magnitude. 

Our initial model extends up to the photosphere. The 
photospheric region is characterised by the ionisation re- 
gions of Hell and Hel/H, with "bumps" in the opacity pro- 
file. In addition, the surface region is characterised by the 
lowest values of the pressure scale height within the whole 
stellar structure. A non-uniform radial mesh would be nec- 
essary to solve the different scales present in the star. Our 
ID stellar structure code uses such an adapted radial grid, 
which ensures a very good resolution of all features through- 
out the star. This is illustrated in the left panel of Fig. [SJ 
where we have sampled the initial ID model to create a 
radial mesh with a suitable number of grid points (200 in 
this case). 

However, when considering multi-dimensional calcula- 
tions, this grid presents a flaw: extension in the tangential 
direction yields intractable aspect ratio of the cells. This 
is illustrated in the right panel of Fig. |8j which shows the 
radial profile of the cell aspect ratio when the former ra- 
dial mesh is simply extended in the tangential direction by 
putting 256 cells over a tangential angular width of 7r/2. 
Inspection of Fig. [5] shows that a constant number of cells 
per scale height translates into a large variation of the cell 



aspect ratio. In the particularly demanding region near the 
photosphere the aspect ratio drops below 10~ 2 , which is 
unacceptable for performing numerical simulations. 

It is therefore obvious that to perform multi- 
dimensional simulations from the deep stellar interior to 
the surface, one has to make a compromise on the resolu- 
tion when using a single grid. A solution to this problem 
would be to use a system of nested grids where the grids are 
refined in both the radial and tangential directions as mov- 
ing near the surface of the star. Here, we choose to consider 
a computational grid with an uniform grid spacing in both 
Ar and A9 and apply a special treatment for the surface 
layers (see below). Due to spherical coordinates, the aspect 
ratio of the cells decreases with the radius. The radial step 
is chosen so that the cells in the middle of the domain are 
square. The right panel of Fig. [5] shows that our cells have 
an aspect ratio of 3 at the bottom of our domain and 0.6 at 
the top of the domain. As a result, the surface of the star 
is not resolved well, as shown by the number of grid points 
per pressure scale height (see left panel of Fig.[H]). 

5.2.2. Simulation setup 

Since we cannot accurately describe the surface of the star, 
we choose to apply a heating function that drives an isother- 
mal region in the outer, badly resolved, region of the star. 
To do so, we include in our energy equation a heating source 
term on the right-hand side: 



-pc v (T - T )/(r)/r h eati] 



(42) 



where /(r) is a spatially varying function that is equal 
to 1 above a given radius r c and zero elsewhere, with a 
smooth transition in between, locating is the typical cooling 
time scale and To is the forcing temperature. This method 
has been used for th e "star-in-box" simulation presented in 
Dobl er et al.l (|2006l ) (see also references therein). The pa- 
rameter To is chosen to be the temperature of the gas at r c 
in the initial model of the star. For the model presented in 
this section, r c — 0.95i?*, To = 32 750 K, and Treating = 10 4 
s. For comparison, the value of the radiative time scale (as 
defined in Eq. l31[) decreases exponentially from ~ 3.5 x 10 5 
s to ~ 3.5 x 10 3 s in this region so our value of Theating 
insures that the source term acts on a short enough time 
scale to keep the outer region in a roughly isothermal state. 

For the value of To used in our model, we do not have the 
Hel/H ionisation region near the surface. Also we prevent 
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Fig. 10. Evolution of the total kinetic energy (top panel) 
and CFL numbers for the giant star run (lower panel, solid 
red: radiative; dashed blue: hydro; dotted green: advection). 
The vertical dashed lines show the time interval on which 
the averaged fluxes shown in Fig. 1121 were computed. 
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Fig. 12. Radial profiles of averaged fluxes for the giant 
star run (dashed-dotted line: radiative flux; dashed line: 
enthalpy flux; dotted line: kinetic energy flux; solid line: 
total flux). 



any convection motion in this region since the isothermal 
structure is stably stratified. There is essentially no radia- 
tive flux throughout this region, because the energy enter- 
ing this region is evacuated by the source term instead. 

This treatment of the surface layers is very basic, and is 
only used as a preliminary solution to produce our first re- 
sults. We are currently investigating other more physically 
based methods to treat the surface layers. 

5.2.3. Results 

We now present the results of the run (summarised in 
Table 0} . Figure [TU] shows the evolution of the kinetic en- 
ergy and CFL numbers. As in the previous section, we see 
that, during the initial phase of relaxation/development 
of convection, we reach high CFL numbers, with the hy- 
dro CFL number reaching a maximum value around 10 3 . 
As the convective instability develops, we see a quick and 
strong rise in the kinetic energy in the domain, followed 
by a slower decrease as the model relaxes toward equi- 
librium. After 500 days, the kinetic energy tends toward 
a constant value, indicating that a quasi steady-state is 
reached. The time step remains roughly constant with a 
mean value of At = 2.5 x 10 4 ± 5 x 10 3 s, corresponding to 
CFL hydro = 100 ± 20, and CFL rad = 7.6 ± 1.5. The second 
value indicated corresponds to the standard mean devia- 
tion. Therefore in this kind of model, we see that the ra- 
diative time scale is much lower than the hydro time scale. 
Figure [TT] shows a snapshot of the convective motions in 
the star. 

Wc computed the averaged fluxes defined in Sect. 15.11 
between t = 3000 days and t — 4000 days (corresponding 
to 3452 time steps) . The resulting radial profiles are shown 
in Fig. [T3] We see that thermal equilibrium has not been 
reached yet, since the total flux throughout the envelope 
is not constant. We see a significant upward enthalpy flux 
in the ~ 40% outer radius of the star. Also, we see a very 
important downward kinetic energy flux in the convective 
zone. Around r ~ 0.5i?*, we observe strong bumps both 
in the radiative and in the enthalpy flux. This region cor- 



10 s 




10 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 
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Fig. 13. Radial profiles of the CFL time steps for the giant 
star run (see text). From bottom to top: hydro (dashed 
blue), advection (dotted green), and radiation (solid red). 
The horizontal dash-dotted line indicates the time step of 
the snapshot under consideration. 

responds to the interface between the convective zone and 
the radiative zone. When downdrafts cross the interface, the 
low efficiency of radiation prevent them from warming up, 
and as a result, they are subject to a very strong braking by 
buoyancy forces, which prevent them to significantly over- 
shoot in the radiative zone. This is why the interface region 
is strongly perturbed whereas the inner radiative zone re- 
mains unperturbed. Finally, we checked the range of Mach 
numbers in the domain. The average Mach number of the 
flow extends from ~ 10~ 3 at the bottom of the envelope to 
~ 0.15 at the top of the convective zone. However, Mach 
numbers in the range 0.30 — 0.40 can be reached locally, 
in the downdrafts. We stress that such a range of Mach 
number cannot be addressed by anelastic codes, which are 
restricted to low Mach number flows (M s < 10 -2 — 10 -1 ). 

Figure [T3] shows the radial profiles of the different CFL 
time steps (as in Fig. [5]). As already mentioned above, we 
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Table 4. Numerical parameters for the cold giant run (same as Table [3]). 



Resolution Domain 


TKH At hydro 


A.Ct'L 


Final time 


Time steps Wall Time 


Mem. 


217 x 256 [0.2,1.] x [7r/4,37r/4] 


2200 yr 250 s 


3.2 x 10 a s 


4235 d 


16216 8 w 


2.2 Gb 



t= 4.310e + 03 days 




-0.6 -0.4 -0.2 0.0 0.2 0.4 0.6 



Fig. 11. Snapshot of the cold giant star simulation. The axis are Cartesian coordinates normalised by the stellar radius 
and with the origin at the star centre. 



see that the radiative CFL time step is much longer than 
the hydro time step; however, the strong decrease in the ra- 
diative time scale in the isothermal surface layer comes from 
the decrease in density in this region. We can expect that, 
if the surface regions were better resolved, the radiative 
time scale would drop to a lower value than the hydro time 
step, and we would recover locally a situation similar to 
the A-type star model. Here again the shortest hydro CFL 
time step is found at the bottom of the envelope, therefore 
we reach a higher hydro CFL number than for the A-type 
star model because our model extends deeper in the star. 
Finally, as discussed in the previous section, we see that 
the actual time step is roughly equal to the lowest value of 
the advective time scale, which is located in the convective 
zone. In this case, since the envelope extends deeper, our 
implicit scheme allows for reaching high values of the hydro 
CFL number by removing the limitation coming from the 
deeper layers of the star. 

The results presented in this section show that it is 
meaningful to adopt a fully implicit approach for comput- 
ing models of deep convective zones. Indeed, we show here 
an example of a situation where substantial CFL numbers 
could be reached by overcoming the time step limitation 
imposed by deep layers of the stellar envelope. This is pos- 
sible since the flow velocity is low in those regions. In the 
surface convective region, however, we are still prone to the 
time step limitation imposed by the flow movement across 
grid cells. 



Finally, the computation of such deep stellar envelopes 
also poses the problem of how to achieve the thermal re- 
laxation of the convective envelope on an acceptable time 
scale. Indeed, the thermal time scale of our numerical do- 
main is 

tkh = ^f- ~ 2200 years, (43) 

whereas the simulation only ran for 4235 days ~ 11.6 years 
of stellar time. This is, however, an unavoidable difficulty 
that is common to all numerical approaches. 

6. Conclusion 

We have presented a new numerical code devoted to stellar 
interiors. The code solves the full set of the compressible 
Euler equations and radiative diffusion with a spatially and 
temporally accurate implicit scheme. 

As a first step, we computed stellar convection following 
two approaches: a local approach for a relatively hot star 
(see Sect. I5.1|) and a global approach for a cold giant star 
(see Sect. 15. 2[) . In the first case, we were able to simulate 
the convective motions in the outer layers of a A-type star. 
For this type of model, we achieve a mean CFL number of 
10 for hydro and of 400 for radiation. The rather low value 
of the hydro CFL number suggests that, in this particular 
case, a mixed approach (explicit hydro -I- implicit radiation) 
would result in a substantial speeding up of the calculation 
compared to a fully explicit or implicit approach. In the 
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second case, we modelled the turbulent convective motions 
that transport energy in the bulk of a cold giant star. Our 
domain includes both the outer convective region and the 
bottom, stable, radiative region. For this model, we could 
reach a mean CFL number of 100 for hydro and of 7 for 
radiation. Such values of the hydro CFL numbers are en- 
couraging and suggest that we might gain computational 
time thanks to our implicit approach. 

The aim of this paper is to show the feasibility of us- 
ing a fully implicit method for the accurate computation of 
time-dependent flows. We stress that our goal is, however, 
not to demonstrate that the code is successful in the sense 
of efficient production. Implicit methods are computation- 
ally much more expensive than explicit methods because 
they need to solve a set of non-linear equations. The effi- 
ciency of an implicit code is essentially set by the efficiency 
of the method used to solve this difficult task. The general 
framework for solving the non-linear system is the Newton- 
Raphson method. For tests and development purposes, we 
chose here a very simple approach where the Jacobian ma- 
trix is computed at each iteration, and where the linear 
system is solved exactly inside the Newton-Raphson loop. 
This results in a robust, but expensive non linear solver. 
A more efficient class of non linear solvers include inexact 
Newton- Krylov methods, which are Newton-Raphson meth- 
ods where an iter ative method is used to solve the linear 
system, see e.g. see lKnoll fc Kev cs 2004 for a review. These 
methods are Jacobian free methods and gain computational 
time by using an iterative method to approximately solve 
the linear system during the Newton-Raphson loop (hence 
the name "inexact"). These methods are very popular in 
computational physics and we plan to implement and test 
them in a future version of the code. Once the most efficient 
non-linear solvers currently available are implemented, we 
will be able to benchmark our code against explicit ones. 

Concerning our numerical method, we found that the 
maximum time step that we can use without being faced 
with convergence issues in the Newton-Raphson procedure 
is closely related to the fluid velocity in the numerical do- 
main. We introduced the advective time scale 



and showed that for the computation of a time-dependent 
flow the time step cannot grow much above the minimum 
value of this time step. The definition of Eq. [44] should be 
adapted in the multi-dimensional case. In other words, our 
numerical method cannot handle fluid movements across 
several cells, which is intuitively and physically understand- 
able. We stress that such a criteria is also suggested by 
accuracy concerns. This limitation on the time step stems 
from the fact that in order to advance the solution over one 
time step, one has to solve a system of non-linear equations 
with a Newton-Raphson procedure. For time step larger 
than the advective time scale, the fluid moves across sev- 
eral cells and the solution changes significantly between two 
time steps. This change is essentially the limiting factor for 
increasing the time step, since the Newton-Raphson proce- 
dure encounters problems when relative changes between 
two time steps are too large. We stress that one could im- 
prove this by finding a better initial guess for the Newton- 
Raphson procedure, instead of simply taking the solution 
of the last time step. In ID implicit simulations, a time 



extrapolation from preceding time steps has been used suc- 
cessfully to produce an initial guess for the Newton proce- 
dure. However, in multi-dimensional flows, this procedure 
does not necessarily gives a satisfactory guess of the solu- 
tion. Indeed our experience did not show any improvement 
by using such an extrapolated guess. 

Our first developments show that using a fully implicit 
code to describe stellar interiors is realistic, providing a 
powerful tool for studying a full stellar structure charac- 
terised by low-to-moderate Mach numbers. Our test per- 
formed on the envelope of a giant star indeed highlights 
the possibility to extend the spatial domain to a large frac- 
tion of the star (80% in radius), where Mach numbers cover 
values from less than 1% to a few tens of percent. More gen- 
erally, this implicit tool will be efficient for computing re- 
laxation toward steady state solutions and very low Mach 
number flows, and to avoid the CFL constraint due to a 
particular location of the numerical domain. In the later 
case, it is possible to have a moderate Mach number flow 
in another region of the domain (this region would then 
be responsible for the limitation of the time step), as illus- 
trated by the red giant case mentioned above. 

Our first developments also stress the numerical diffi- 
culty of modelling a full star with a grid-based code. Since 
the typical length scales (e.g. pressure scale height) vary 
significantly between the deep interior of the star and the 
surface, it is not possible to resolve both the interior and 
the surface layers on a single grid. We plan to overcome 
this difficulty by implementing a system of nested grids re- 
fined in both the radial and tangential directions. Among 
other future developments of the code, we plan to extend 
it to 3D and to implement iterative solvers, rather than a 
direct method, involving parallelisation. Among our first as- 
trophysical applications, we quote (i) the very early stages 
of evolution of low-mass stars, bridging the gap between 
the hydrodynamical phase of star formation and the quasi- 
static pre-main sequence evolutionary phase, (ii) the study 
of convection in pulsating stars, and (iii) the effect of fast 
rotation on convection (in low-mass stars, massive stars). 
But clearly, a wide range of astrophysical problems can be 
studied with such a tool. Based on our experience, we hope 
to achieve these future developments and the first astro- 
physical applications within a reasonable period. We also 
hope that our efforts demonstrate the feasibility of apply- 
ing an implicit multi-dimensional hydrodynamical code to 
stellar interiors and will motivate the computational astro- 
physics community to develop such a tool. 
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Appendix A: Tests cases 

A.l. Linear test cases 
A. 1.1. Advection equation 

In this section we test our numerical scheme with the ID 
scalar advection equation: 



dq dq 
dt dx 



0. 



(A.l) 



(A.2) 



We use a smooth profile as an initial condition: 

f(x) = sin(a;), x € [0, 2ir]. 

The advection velocity a is set to 1. We advect the sin 
wave for one unit of time and the resulting numerical so- 
lution qi + i/2 is compared with the exact solution q® +1 / 2 = 
f{x i+ i/2 — 1) by computing the L\ and errors: 



,.o 



\q-q 



max|<? i+1/2 - <? i+1 . 2 | 



(A.3) 
(A.4) 



For a given time step and mesh size, the CFL number 
is defined as aAt/ Ax. We use an uniform discretisation in 
space and perform a sequence of runs with different res- 
olutions ranging from 50 to 800 grid points by successive 
increase by a factor of 2. This is done for three values of 
the time step in order to investigate a wide range of CFL: 
At = 10 _1 , 10~ 2 , 10~ 3 . We use the Van Leer limiter and 
the Crank-Nicholson method, which are both second-order 
methods. 

Such a grid of simulations allows to study the conver- 
gence properties of our code. The actual order r of the 
scheme between two consecutive resolutions is computed 
as 



log 



\q-q 



011(1) 



\q - q 



011(2) 



/log 2, 



(A.5) 



where (2) corresponds to a resolution twice higher than (1). 
The results are shown in Table IA.11 IA.21 and IA.31 

Before discussing the results, it is useful to derive a re- 
sult that will help for their interpretation. Going back to 
Eq. (|A.1I) , the simplest implicit scheme that can be designed 
combines the first order backward Euler method with the 
first-order upwind method: 



q, 



n+l 



aAt 
Ax 



n+l 

qt-i 



(A.6) 



For this scheme, it is straightfo rward to deri ve the so- 
called "modified equation" (see e.g. lHirschlll990l ): 
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Table A.l. Linear advection, At = 10 _1 . 



N 


Li-error 


Rate 


Loo-error 


Rate 


CFL 


Wall Time (s) 


50 


1.380e-02 




1.152e-02 




7.80e-01 


0.9 


100 


3.457e-03 


2.00 


2.587e-03 


2.16 


1.58e+00 


1.8 


200 


3.138e-03 


0.14 


7.491e-04 


1.79 


3.17e+00 


3.5 


400 


3.252e-03 


-0.05 


8.115e-04 


-0.12 


6.35e+00 


7.7 


800 


3.308e-03 


-0.02 


8.269e-04 


-0.03 


1.27e+01 


27.5 



Table A. 2. Linear advection, At — 10 2 . 



N 


Li-error 


Rate 


Loo-error 


Rate 


CFL 


Wall Time (s) 


50 


1.890e-02 




1.292e-02 




7.80e-02 


3.8 


100 


4.481e-03 


2.08 


4.882e-03 


1.40 


1.58e-01 


7.5 


200 


1.027e-03 


2.13 


1.830e-03 


1.42 


3.17e-01 


15.3 


400 


2.094e-04 


2.29 


6.732e-04 


1.44 


6.35e-01 


34.5 


800 


3.996e-05 


2.39 


1.936e-04 


1.80 


1.27e+00 


134.2 



Table A. 3. Linear advection, At = 10 



N 


Li-error 


Rate 


Loo-error 


Rate 


CFL 


Wall Time (s) 


50 


1.895e-02 




1.293e-02 




7.80e-03 


35.2 


100 


4.530e-03 


2.06 


4.892e-03 


1.40 


1.58e-02 


68.5 


200 


1.074e-03 


2.08 


1.845e-03 


1.41 


3.17e-02 


134.4 


400 


2.569e-04 


2.06 


6.918e-04 


1.42 


6.35e-02 


273.9 


800 


6.102e-05 


2.07 


2.569e-04 


1.43 


1.27e-01 


1016.7 



which shows that formula I A. 6 1 solves an advection-diffusion 
problem since the leading term of the truncature error is 
physically equivalent to a viscous term with a numerical 
viscosity coefficient D = |-(Ax + At). Usually, with explicit 
schemes in mind, it is supposed that the ratio Ax/ At is kept 
constant (i.e. constant CFL number) so that D oc Ax or 
D on At, which is equivalent to the statement that scheme 
IA.6I is first order in time and in space. On the other hand, 
when an implicit scheme is used, we cannot make such an 
assumption for the ratio Ax/ At, as the CFL number is 
in principle arbitrary. In our discussion we use the expres- 
sion of this numerical viscosity coefficient to interpret our 
results. 

It appears that for CFL values above 1 (e.g. At = 10 _1 , 
see Tab. IA. 1|) . the numerical error remains roughly con- 
stant when the resolution is increased. This suggests that, 
in this case, the time truncation error dominates the spa- 
tial one. With the expression of the numerical viscosity de- 
rived above, we are in a regime where D = #At = cte, so 
when the resolution is increased the numerical dissipation 
does not decrease. Another way of interpreting this result is 
that the solution is advected over more than one cell width 
during one time step so that a strong numerical damping 
results. On the other hand, for low values of the CFL num- 
ber (e.g. At = 10~ 3 , see Table |A~3|) . the spatial truncation 
error dominates, and we recover the second-order conver- 
gence of our spatial scheme in the Li norm. We do not 
see an accumulation effect of the round-off errors for very 
small time step. Again using the numerical dissipation de- 
rived above, one has for a first-order scheme D = §Ax. 
We expect that for our second-order scheme, a similar cal- 
culation would lead to D oc Ax 2 , in agreement with the 
observed tendency. A lower rate of convergence measured 



in the Loo is characteristic of the action of the limiter on 
the smooth extrema of the solution. This is a well known 
drawback of TVD schemes: their order decreases also at 
smooth extrema of the solution. Finally for At = 10~ 2 , we 
are in the regime where Ax ~ At, so when the mesh size 
is decreased, the numerical dissipation decreases in a non 
trivial fashion. For instance, considering the first-order nu- 
merical dissipation derived above, decreasing the mesh size 
by a factor of two would result in a decrease in the numer- 
ical dissipation by a factor of fjp^ , where A = Ax/ At. As 
a result the scheme does not behave as having a specific 
order in this case. 

To summarise these results, we stress that the error 
in the numerical solution of partial derivative equations 
depends on both the spatial and temporal discretisation. 
When the temporal (resp. the spatial) truncation errors 
dominates, i.e. for very large (resp. very small) CFL num- 
bers, decreasing the spatial (resp. temporal) step does not 
lead to more accurate results. It is therefore essential to 
choose both the mesh size and the time step in order to op- 
timise the accuracy and the computational cost. Finally, we 
stress that the advection equation is not a suitable prob- 
lem for implicit methods, since the stability criterion for 
an explicit scheme, namely a At < Ax, is also an accu- 
racy criterion. The optimum accuracy /computational cost 
is therefore obtained for CFL number of order unity. Any 
advection computation at high CFL numbers suffers from 
serious numerical damping and lack of accuracy. 

A. 1.2. Diffusion equation 

In this section, we test our diffusion scheme with the simple 
linear diffusion equation: 



17 



M. Viallet et al.: Implicit calculations of stellar interiors 
Table A.4. Linear diffusion, At = 1(T 2 . 



N 


Li-error 


Rate 


ioo-error 


Rate 


CFL 


Wall Time (s) 


50 


2.072e-03 




1.824e-03 




1.50e+00 


0.4 


100 


5.237e-04 


1.98 


3.659e-04 


2.32 


6.13e+00 


0.6 


200 


8.663e-04 


-0.73 


8.220e-04 


-1.17 


2.48e+01 


1.0 


400 


9.727e-04 


-0.17 


9.487e-04 


-0.21 


9.95e+01 


2.0 


800 


9.999e-04 


-0.04 


9.802e-04 


-0.05 


3.99e+02 


6.6 



Table A.5. Linear diffusion, At = 10" 3 . 



N 


Li-error 


Rate 


Loo-error 


Rate 


CFL 


Wall Time (s) 


50 


2.825e-03 




2.812e-03 




1.50e-01 


2.6 


100 


6.874e-04 


2.04 


6.743e-04 


2.06 


6.13e-01 


4.9 


200 


1.629e-04 


2.08 


1.586e-04 


2.09 


2.48e+00 


9.2 


400 


3.419e-05 


2.25 


3.158e-05 


2.33 


9.95e+00 


18.3 


800 


4.973e-06 


2.78 


3.406e-06 


3.21 


3.99e+01 


64.4 



Table A.6. Linear diffusion, At = 10" 4 . 



N 


Li-error 


Rate 


Loo-error 


Rate 


CFL 


Wall Time (s) 


50 


2.833e-03 




2.823e-03 




1.50e-02 


24.8 


100 


6.963e-04 


2.02 


6.847e-04 


2.04 


6.13e-02 


47.9 


200 


1.717e-04 


2.02 


1.690e-04 


2.02 


2.48e-01 


91.4 


400 


4.265e-05 


2.01 


4.193e-05 


2.01 


9.95e-01 


183.7 


800 


1.059e-05 


2.01 


1.038e-05 


2.01 


3.99e+00 


640.5 



For an initial Dirac distribution of amplitude Q at 
x=0, the analytic solut ion of Eq. (|A.8p reads (see e.g. 
iMihalas fc Mihalas|[l999T) as 

T(z,t) = ^L=exp(-^-). (A.9) 

The problem is solved numerically on a uniform grid span- 
ning the domain [—2,2]. For a given mesh size and time 
step, the CFL number is defined as x At/ Arc 2 . We use our 
second-order spatial scheme for diffusion, together with the 
second-order Crank-Nicholson scheme for time stepping. 

For simplicity, we consider here x = Q = 1. Since it 
is not possible to represent a Dirac distribution exactly on 
a discrete mesh, we adopt as an initial condition the solu- 
tion of Eq. (|A.9[) at t=0.025. The solution is evolved un- 
til t=l, and we compare the numerical solution with the 
analytic solution, as done in the previous section for the 
linear advection problem. Again we perform a sequence of 
grid resolutions ranging from 50 to 800 grid points. This is 
done for three different time steps At = 10~ 2 , 10~ 3 , and 
10~ 4 , yielding a wide range of CFL numbers. The results 
are shown in Tables IA.41 IA.5I and IA.6I 

As in the previous section, when the time step is large 
(At = 10~ 2 , see Table [A~4|) . the temporal truncation error 
dominates, and the numerical error does not decrease when 
the resolution is increased. On the other hand, for small 
time steps (e.g. At = 10~ 4 , see Table IA.4|) . the spatial 
truncation error dominates, and we recover a second-order 
convergence. In this case, both the L\ and the L m show 
second-order behaviour, since there is no limiting process 
in the construction of the interface fluxes. Again in the 
transition regime (At = 10~ 3 , see Table |AT4|) . we do not 
observe a well-defined order. Such behaviour is likely for a 
similar reason to the advection case. 

Here, an implicit scheme is more efficient than an ex- 
plicit scheme since, for a given resolution, for instance 
N=800, the run with At = 10~ 2 has a solution that is 



more accurate (by more than a factor of 10 in both norms) 
than the run with At = 10~ 3 (closer to the explicit time 
step), and has a computation time that is ten times shorter. 
In this case, the optimum ratio between accuracy and com- 
putational cost is obtained for CFL numbers larger than 1, 
because here the CFL condition is not an accuracy crite- 
rion. 

A.2. Non-linear test cases 
A.2.1. Sod test 

A typical test for non-linear hydro is the ID Sod shock-tube 
test. This is a Riemann problem with the following initial 
right and left states: ph = 1, Pl = L PR = 0.125, pr = 0.1; 
the gas is initially at rest. The solution of the Sod problem 
is computed until t = 0.25 on the domain [—0.5, 0.5], where 
the interface between the initial states is located at x = 0. 
The number of grid points is 400. For simplicity and since 
the goal is to show how our code deals with shocks, we set 
the time step to keep the hydrodynamic CFL equal to 1. 

The left panel of Fig. lA.ll shows the result of the compu- 
tation. One can easily identify the shock front, the contact 
discontinuity, and the rarefaction wave. The location and 
amplitude of these features are in excellent agreement with 
the analytical solution. However, one can see that spuri- 
ous numerical effects affect the solution. Oscillations are 
produced at the location of the shock, despite the use of a 
limiter. An over/under shoot is also visible at the tail of the 
rarefaction wave. Solving the total energy equation and/or 
limiting the primitive variables (i.e. pressure) could improve 
these results. Otherwise, a possible cure is to introduce an 
artificial viscosity of the form: 

z/ = CAxc 2 , (A.10) 

where c s is the sound speed, Ax the mesh size, and C a 
dimensionless parameter that has to be fine-tuned to obtain 
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acceptable results. The right panel of Fig. IA.1I shows the 
result of the computation when C = 0.5. We stress that 
we do not encounter shocks in our current calculations of 
stellar models. 

A. 2. 2. Barenblatt test 

To test our diffusion scheme in the non-linear regime, we 
performed various computations of the propagation of non- 
linear conduction wave, also known as the Barenblatt prob- 
lem. We solve the equation 



(A.ll) 



dT _ 

dt dx\ dxJ' 
which is a diffusion equation with a non-linear diffusivity 
law x — ctT@ . 

As in the linear diffusion problem, we assume that an 
instantaneous pulse of energy is released at x = 0. Let Q 
be the amplitude o f the temperature pulse re leased in the 
medium. Following Mihalas & Mihalas (1999), the analyti- 
cal solution of the problem is 



with 

/(0 



208 + 2) 
= 0, otherwise. 

£o is defined by 



(A.12) 



(-^-) 

\2(B + 2)J 



VP r 



for£<£o (A.13) 
(A.14) 



(A.15) 



- (/3 + 2)^+ 1 2 1 -^r(^ + l//3) i V0+2 

. /3tt^/ 2 r(l//3) - 

The problem can be cast in dimensionless variables, defined 
in the following way: 



t = aQ p t 
x = x/£ a 

f = t/q/(o). 

In these variables, Eq. (IA.11[) becomes 

8 d (-odT 



dT 



2(8 + 2)dx\ dxJ' 
and Eq. (|A.12|) now reads as 



T(x,t) = t 



-1/0+2 



1 



(fL/0+2) 



1/0 



(A.16) 
(A.17) 
(A.18) 



(A.19) 



(A.20) 



The advantage is that, in the dimensionless form, the 
only parameter left is 8. We performed runs for 8 = 1, 3, 
5,7, going from mildly to strongly non-linear regimes. The 
problem is solved in the dimensionless form for the domain 
[-1.5,1.5]. Our spatial discretisation uses 100 grid points, 
and again the time step is set so that the CFL number is 
equal to 1. As we cannot represent a 5 pulse on a discrete 
grid accurately, the initial condition in each case is taken 
to be the analytical solution (IA.20j) at t = 0.1. Figure lA~2l 
shows the results at different times and demonstrates the 
ability of our code to resolve and track the heat front even 
for strongly non-linear cases. 



Appendix B: Formulae 

We recall the geometrical factors for the scalar, cell centred, 
quantities: 



2tt 



i+1 



■)(cos( 



COS( 



3+1 j 



S lj+i/2 = 2nrj (cos 6j - cos^+i) 
^ 1/2j =™^(r 2 +1 -r 2 ). (B.14) 

Geometrical factors for the radial velocity 

27T 

Kj+1/2 = y(ii/ 2 ~ r^ 1/2 ){cos9j - cos^+i) 

"ST-1/2J+1/2 = 2«ii/2(a» fl 3 - COS^+l) 

Sf J = n&i6jtf +1/t -j_ 1/2 ). (B.15) 

Geometrical factors for the tangential velocity 
2ir 

Vi + l/2,j = y K+l ~ r i)( COse j-l/2 - C0s6» i+1/2 ) 

= 27rr I 2 (cos6» J _ 1/2 - cos 6> J+1/2 ) 

Sf+i/2,j-i/2 = 77 sin ^-1/2(^+1 - r\). (B.16) 

The complete set of ordinary differential equations resulting 
from our spatial discretisation in spherical coordinates are 
summarised in Tab. IB. II 
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Fig. A.l. Results of the Sod test at t = 0.25 without artificial viscosity (left panel) and with C — 0.5 (right panel). The 
dots are the numerical solution, the solid line is the analytical solution. The time step corresponds to CFLi iy d ro = 1- 




Fig. A. 2. Results of the Barenblatt problem for /3 = 1, 3, 5, 7. The solution is shown at i = 0.1, 1, 2 for (3 = 1 and at 
i = 0.1, 1, 5 for (3 — 3, 5, 7, from upper to lower curves. The crosses are the numerical solutions, the solid lines are the 
analytical solutions given by Eq. (|A.20[) . 
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^Pi+l/2,j+l/2V j+ l/ 2 j + l/2 - -(ST + ij+i/ 2 Pi+lj+l/2Ui+lj+l/2 - Si,j + l/2Pi,j + l/2K,j + l/2 

+'-'i+l/2j+lPi+l/2j+l'" j+ i/2j+l — &i Pi+l/2,jUi+l/2,j) 
Pi+l/2,j + l/2e i+ l/2,j + l/2 1 K + i/2,j + l/2 = -(-ST+lj+l/2Pe i+ljJ - +1/2 «i+l,j+l/2 ~ 'ST,j+l/2Pe jjJ - +1/2 Mi,j+l/2 



+'S'i+l/2j+lPe i+ i/ 2j j + i«i+l/2,j+l — Si+l/2,jP e i+l/2,j u i+l/2,j) 



~Pi+l/2,j+l/2{S i+ i t j + i/2 u i+l,j+l/2 ~ •S'i,j+l/2' lt i,j + l/2 + 'S'i+i/2, :) '+lttj +1/ '2,j + l — >S'i+l/2,j' lt i+l/2,j) 



i+lj+l/2 = 



+3/2,j+l/2 



i+l/2,j + l/2 



P K i+l j+1/2 



t; 



i+l/2j+l = 



A i+ ir 

+ l/2,j+3/2 _ r /+l/2,j + l/2 



+l/2,j+l/2 



-l/2j + l/2 



P K i,j+l/2 



Air 



P K i+l/2,j+l 
d_ 

dt 



r i+1/2 A j+1 e 

Pi,j+l/2Ui,j+l/2Vi,j+l/2 



d _ 9 T > 

P»+l/2,jM»+l/2,j ^+1/2J 



t; 



+1/2 j : 



+ 1/2J + 1/2 i i+l/2,j-l/2 ' 



P K i+i/2,j r i+1/2 AjG J 

{ S i+l/2,j + l/2pU r i+1 / 2 ,j + l/2 u i+l/2,j + l/2 
- ST_i/ 2 ,j + l/2p« r 'i-l/2,j+l/2^i-l/2,j+l/2 

+ S , i, J - + ipW ij - +1 iii,j + i - SijfWijUij) 

Pi+l/2,j+l/2 ~ Pi-l/2,j+l/2 ,J ,y 

T— ^,3 + 1/2 + Pi,j+l/29i,j+l/2 Vi, 3 + l/2 

Pi,] + l/2{u e ij + l/2 f 



n 



-(SI 



a — r r-iT a — r 

+1 , 3 pu v i+1:j u i+lij - •S,. ( pu",. ( ii i 



+ Sf 



+1/2J+1/2PU i+l/2,j+l/2 u i+l/2,j+l/2 
+ l/2,i-l/2p" 9 l+ l/2,j-l/2'"i+l/2,j-l/2) 



f<+l/2j+l/2 ~ Pi+l/2,j-l/2 .-. 



(B.l) 
(B.2) 



(B.3) 



(B.4) 



r i+ i /2 Aj0 



^+1/2 j - ^»+l/2,j 



-r 8 
Pi+l/2,jUi+l/2,j u i+l/2,. 



r i+l/2 



(B.5) 



P K i,j+l/2 
1 



Pi+l/2,j 

i+l/2,j+l/2 

-r 

«I+l/2j 



1.1 1 

o( + 

^ Pi-l/2,j+l/2«i-l/2j+l/2 Pi+l/2j+l/2«»+l/2j+l/2 

1 1 1 

o v 



P K i+l/2,j 2 Pi+i/2j-l/2«i+l/2,j-l/2 Pi-l/2,j+l/2K»+l/2j+l/2 

Pij+1/2 = (Pi-l/2j+l/2Vi-l/2j+l/2 + Pi+l/2,j + l/2V r i + l/2,j + l/2)/V r l , J + l/2 
(Pi+l/2,j-l/2 , K + l/2 i j + l/2 + Pi+l/2,j-l/2Vi + l/2,j-l/2)/V i+ l/2j 



2( W L+l/2 + «i+l,j+l/2) 
2( U i,j-l/2 +WL+I/2) 



(B.6) 

(B.7) 

(B.8) 
(B.9) 

«L = ^(Ui-l/2,3 + u 1+l/2,j) (B.10) 
«i+l/2,j+l/2 = ^(*4+l/2J +«i+l/2,j+l) (B.ll) 



ij + 1/2 



, -6 



3 + 1) 



(B.12) 
(B.13) 



Table B.l. Complete system of ordinary differential equations resulting from our finite volume discretisation on a 
staggered grid in spherical coordinates. 
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